rm(list=ls())
###
### Replication material for Binetti & Steinwand,
### "Aid targeting in post-conflict Nepal",
### published in Conflict Management and Peace Science
### 


#Set your path to the replication ordner 
	#containing the replication data "BinettiSteinwand_CMPS_data.rds"

#setwd(<mypath>)
setwd("/Users/martin/Dropbox/Documents/ProjectAidPostConflict/PCATarget_shared/Data/CMPS_replication")


#Load replication data
load("BinettiSteinwand_CMPS_data.RData")

#create data files
df_nepal <- dat_all$df_nepal
Nepal.map_sf <- dat_all$Nepal.map_sf
df_nepal_dist <- dat_all$df_nepal_dist
dist_inst <- dat_all$dist_inst
nb <- dat_all$nb
lw <- dat_all$lw
lw_bin <- dat_all$lw_bin
KillRatio <- dat_all$KillRatio
barracks <- dat_all$barracks
Nepal_map_agg_sf <- dat_all$Nepal_map_agg_sf
Nepal_map_dist_sf <- dat_all$Nepal_map_dist_sf
crs_raw <- dat_all$crs_raw
M_TotAid_SAC <- dat_all$M_TotAid_SAC
M_SocAid_SAC <- dat_all$M_SocAid_SAC
M_InfrAid_SAC <- dat_all$M_InfrAid_SAC


###Please load the following packages

#library sf should be updated to 1.0 or higher to enable spherical mapping: https://r-spatial.org/r/2020/06/17/s2.html
library(tidyverse)
library(sf)
library(raster)
library(stargazer)
library(geodata)
library(geojsonio)
library(ivreg)
library(spdep)
#for ability to customize stargazer output
#devtools::install_github("ChandlerLutz/starpolishr")
#source: https://rdrr.io/github/ChandlerLutz/starpolishr/
library(starpolishr)
library(texreg)
library(marginaleffects)

library(olsrr)
library(sandwich)
#library(lmtest)
library(car)

# library(devtools)
# # devtools::install_github("chadhazlett/sensemakr")
# library(sensemakr)
library(spatialreg)

library(sphet)
library(plm)

library(tmap)

#library(readxl)


#########
#Figure 1
Distmap <- st_as_sf(Nepal.map_sf)
df_nepal <- as_tibble(df_nepal)

#wanted: total aid, out-of-area killings, in-area killings, vote share Madhesi
#total aid
aid_breaks=c(0.50, 20, 1100)
aid_plot <- ggplot(data = df_nepal) +
    geom_sf(aes(geometry=geometry, fill=aid_all_tot_pc+.01)) +
    coord_sf(label_axes = "----") +
    scale_fill_gradient2(name = "USD per capita,\nlog scale", trans = "log", 
      breaks=aid_breaks, low="green", high="blue") +
    labs(caption = "(a) Total aid committed") +
    theme(plot.caption = element_text(hjust = 0))

ggsave("aid_map.pdf", width=6, height=3, units="in")


aid_social_plot <- ggplot(data = df_nepal) +
    geom_sf(aes(geometry=geometry, fill=aid_socl_tot_pc+.01)) +
    coord_sf(label_axes = "----") +
    scale_fill_gradient2(name = "USD per capita,\nlog scale", trans = "log", 
      breaks=aid_breaks, low="green", high="blue") +
    labs(caption = "(b) Social sector aid committed") +
    theme(plot.caption = element_text(hjust = 0))

ggsave("aid_social_map.pdf", width=6, height=3, units="in")

aid_infra_plot <- ggplot(data = df_nepal) +
    geom_sf(aes(geometry=geometry, fill=aid_infr_tot_pc+.01)) +
    coord_sf(label_axes = "----") +
    scale_fill_gradient2(name = "USD per capita,\nlog scale", trans = "log", 
      breaks=aid_breaks, low="green", high="blue") +
    labs(caption = "(c) Infrastructure aid committed") +
    theme(plot.caption = element_text(hjust = 0))
ggsave("aid_infra_map.pdf", width=6, height=3, units="in")


#########
#Figure 2

#code in Nepal_code_barracks.rmd

# transform the dataset in a spatial object 
DF_battalions_sf <- barracks %>%
	dplyr::select(N_battalions,
		longitude_centroid,
		latitude_centroid) %>%
	filter(N_battalions >0) %>%
	st_as_sf(coords = c("longitude_centroid", "latitude_centroid"),
                             crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 + no_defs")

png("map1.png",
    width=5000,
    height=3286,
    res=600)

###Nepal_map_dist_sf does not exist 

tm_shape(Nepal_map_dist_sf)+ # selct polygon to plot
  tm_borders(lwd=1,         # select thickness of the line  
             fill = "yellowgreen", # fillour of borders
             lty = "solid")+
  tm_shape(Nepal_map_agg_sf)+ # selct polygon to plot
  tm_borders(lwd=2,         # select thickness of the line  
             # fill = "black", # fillour of borders
             lty = "solid")+
  # tm_fill(fill = "#ad9a00",
  #         fill_alpha = 0.5)+ # type of the line of the border
  tm_shape(DF_battalions_sf)+ 
  tm_symbols(fill="#006aad",
             fill_alpha = 1,
             size=0.5)+
  tm_credits("Source: Nepal Conflict Archive", 
             position=c("left", "bottom"),
             bg.fillor = "white",
             bg.alpha = 0.1)+
  tm_add_legend('symbol', 
                fill = c("#006aad"),
                fill_alpha = 1,
                shape = c(21),
                labels = c("District centroid"),
                is.portrait = T,
                title="Location of barracks")
  # tm_layout(bg.fillor = "white",
  #           frame = FALSE,
  #           main.title.position = "center",
  #           main.title.size= 1.2,
  #           legend.title.size = 1, 
  #           legend.text.size = 0.8,
  #           title.fontface="bold",
  #           legend.title.fontface= "bold",
  #           legend.height=10)

dev.off()

#########
#Figure 3

# maps at district level 
## merge in data
df_nepal_dmap <- merge(x = df_nepal_dist, y = KillRatio, by = "district_code")
df_voting <- df_nepal %>%
  group_by(district_code) %>%
  summarize(vote_mao = max(Vote.share.Maoist), 
    vote_madesh = 
    max(Vote.share.Terai_madesh+Vote.share.Madeshi_jana),
    CivGovShare = first(CivGovShare))
df_nepal_dmap <- merge(x = df_nepal_dmap, y =df_voting , by = "district_code")
df_nepal_dmap <- df_nepal_dmap %>%
  st_as_sf %>%
  dplyr::mutate(kill_non_mao=1-kr, totkill = ifelse(is.na(totkill), 0,
    totkill))

civkilled_plot <- ggplot(data = df_nepal_dmap) +
  geom_sf(aes(geometry=geometry, fill=CivGovShare)) +
  coord_sf(label_axes = "----") +
  scale_fill_gradient(name = "Share of\ncivilian deaths",
  low = "light gray", high="dark red") +
  labs(caption = "(a) Civilians killed by government") +
  theme(plot.caption = element_text(hjust = 0))

ggsave("civkilled_map.pdf", width=6, height=3, units="in")


maovote_plot <- ggplot(data= df_nepal_dmap) + 
  geom_sf(aes(geometry=geometry, fill=vote_mao)) +
  coord_sf(label_axes = "----") +
  scale_fill_gradient2(name = "Vote share CPN(M)",
  high="orange") +
  labs(caption = "(b) Voter support for CPN(M)") +
  theme(plot.caption = element_text(hjust = 0))
ggsave("maovoteshare_map.pdf", width=6, height=3, units="in")

madhesh_plot <- ggplot(data = df_nepal_dmap) +
  geom_sf(aes(geometry = geometry, fill = vote_madesh)) +
  coord_sf(label_axes = "----") +
  scale_fill_gradient2(name = "Vote share\nMadhesi parties",
    high = "purple") +
  labs(caption = "(c) Voter support for Madhesi parties") +
  theme(plot.caption = element_text(hjust = 0))
ggsave("madheshvoteshare_map.pdf", width=6, height=3, units="in")



########
#Table 1 Summary of variables in the analysis
sumdat <- df_nepal %>% 
    dplyr::select(aid_all_tot_pc_ln, aid_socl_tot_pc_ln, aid_infr_tot_pc_ln,
    CivGovShare, Vote.share.Maoist, Vote.share.Terai_madesh, Vote.share.Madeshi_jana, totkill, 
    v4composites_calibrated_201709.2006.mean, v4composites_calibrated_201709.1996.mean,
    srtm_slope_500m.none.mean, srtm_elevation_500m.none.mean) %>% 
    mutate(lights_delta=v4composites_calibrated_201709.2006.mean - 
          v4composites_calibrated_201709.1996.mean, 
          v4composites_calibrated_201709.1996.mean=NULL,
          Madeshi_votes = Vote.share.Terai_madesh + Vote.share.Madeshi_jana,
          Vote.share.Terai_madesh = NULL,
          Vote.share.Madeshi_jana = NULL) %>%
    relocate(Madeshi_votes, .after = Vote.share.Maoist) %>%
    relocate(lights_delta, .after = totkill)

star_sum <- stargazer(sumdat, covariate.labels=c("Total Aid per Capita (USD log)", 
    "Social Aid per Capita (USD log)",  "Infrastructure Aid per Capita\\\\(USD log) ",
    "Civilians killed by Govt\\\\\\% of civilian deaths",
    "Maoist Vote Share",
    "Madhesi Parties Vote Share", 
    "Total Fatalities, Average\\\\per Municipality in District", 
    "Change in Economic Activity\\\\1996-2006",
    "Economic Activity 2006", 
    "Average Slope of Area", 
    "Avearge Elevation of Area"),
    title = "Overview of variables",
    label = "tab:vars_summary",
  table.placement = "!htb",
  float = FALSE)

cat(star_sum, file = "summary_alt.tex", sep="\n")

#############
# Tables 2 & 3 OLS and IV results


#Run all models first

##### All aid
#non-IV
F_TotAid_ln_elect_alt <- formula(aid_all_tot_pc_ln ~
  #Political variables
  CivGovShare  + Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean)
      
M_TotAid_ln_elect_alt <- lm(F_TotAid_ln_elect_alt, df_nepal) 

#IV
F_TotAid_ln_elect_iv <- formula(aid_all_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean |
  CivGovShare | N_battalions_50)

M_TotAid_ln_elect_iv <- ivreg(F_TotAid_ln_elect_iv, data=df_nepal) 


######## 
## Social aid

#non-iv
F_SocAid_ln_elect_alt <- formula(aid_socl_tot_pc_ln ~
  #Political variables
  CivGovShare  + Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean)

M_SocAid_ln_elect_alt <- lm(F_SocAid_ln_elect_alt, df_nepal) 


#iv
F_SocAid_ln_elect_iv <- formula(aid_socl_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean |
  CivGovShare | N_battalions_50)


M_SocAid_ln_elect_iv <- ivreg(F_SocAid_ln_elect_iv, data=df_nepal) 


#########
### Infrastructure aid

#non-iv
F_InfrAid_ln_elect_alt <- formula(aid_infr_tot_pc_ln ~
  #Political variables
  CivGovShare  + Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean)

M_InfrAid_ln_elect_alt <- lm(F_InfrAid_ln_elect_alt, df_nepal) 


#iv
F_InfrAid_ln_elect_iv <- formula(aid_infr_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean |
  CivGovShare | N_battalions_50)

M_InfrAid_ln_elect_iv <- ivreg(F_InfrAid_ln_elect_iv, data=df_nepal) 

#compute clustered standard errors
se1_clust <- sqrt(diag(vcovCL(M_TotAid_ln_elect_alt, cluster= ~NAME_2, type = "HC1")))
se2_clust <- sqrt(diag(vcovCL(M_SocAid_ln_elect_alt, cluster= ~NAME_2, type = "HC1")))
se3_clust <- sqrt(diag(vcovCL(M_InfrAid_ln_elect_alt, cluster= ~NAME_2, type = "HC1")))



#for iv models

se1_iv_clust <- sqrt(diag(vcovCL(M_TotAid_ln_elect_iv, cluster= ~NAME_2, type = "HC1")))
se2_iv_clust <- sqrt(diag(vcovCL(M_SocAid_ln_elect_iv, cluster= ~NAME_2, type = "HC1")))
se3_iv_clust <- sqrt(diag(vcovCL(M_InfrAid_ln_elect_iv, cluster= ~NAME_2, type = "HC1")))


#final table non-iv clustered se
table_main_clust <- list(M_TotAid_ln_elect_alt, M_SocAid_ln_elect_alt, M_InfrAid_ln_elect_alt)
star_main_clust <- stargazer(table_main_clust,
  # order = c(1, 7, 8, 2, ),
  dep.var.labels=c("Total Aid", "Social Aid", "Infrastructure\\\\&&&Aid"),
  se = list(se1_clust, se2_clust, se3_clust),
  t.auto = TRUE,
  no.space=TRUE,
  covariate.labels=c("Civilians killed by Govt, \\% of civilian deaths", "CPN(M) vote share",  
    "Madhesi Parties Vote Share", "Total Fatalities", 
    "Change in Economic Activity 1996-2006",
    "Economic activity 2006", "Average Slope of Area", 
    "Average Elevation of Area",
    "Constant"),
  omit.stat=c("ser", "f"),
  title="OLS: Per capita aid by municipality (log USD)",
  label="tab:main_results",
  table.placement = "!htb",
  float = FALSE
#  notes="Robust standard errors. *p$<$0.1; **p$<$0.05; ***p$<$0.001.",
#  notes.append=FALSE,
#  notes.align = "r"
)

star_main_final_clust <- star_insert_row(star_main_clust, "\\textbf{\\textit{Political Support}}\\\\", insert.after = 11)
star_main_final_clust <- star_insert_row(star_main_final_clust, "\\textbf{\\textit{Controls: Reconstruction Need}}\\\\", insert.after = 18)
star_main_final_clust <- star_insert_row(star_main_final_clust, "\\textbf{\\textit{Controls: Economic Need}}\\\\", insert.after = 23)


cat(star_main_final_clust, file = "main_aid_clust.tex", sep="\n")


#final table iv, clustered se
table_iv_clust <- list(M_TotAid_ln_elect_iv, M_SocAid_ln_elect_iv, M_InfrAid_ln_elect_iv)

  #collect diagnostic statistics
summ_tot <- summary(M_TotAid_ln_elect_iv, vcov. = function(x) 
  vcovCL(x, cluster =df_nepal[-M_TotAid_ln_elect_iv$na.action,]$NAME_2, type="HC1"), diagnostics=T)
summ_soc <- summary(M_SocAid_ln_elect_iv, vcov. = function(x) 
  vcovCL(x, cluster = df_nepal[-M_SocAid_ln_elect_iv$na.action,]$NAME_2, type="HC1"), diagnostics=T)
summ_infr <- summary(M_InfrAid_ln_elect_iv, vcov. = function(x) 
  vcovCL(x, cluster = df_nepal[-M_InfrAid_ln_elect_iv$na.action,]$NAME_2, type="HC1"), diagnostics=T)

star_iv_clust <- stargazer(table_iv_clust,
  # order = c(1, 7, 8, 2, ),
  dep.var.labels=c("Total Aid", "Social Aid", "Infrastructure\\\\&&&Aid"),
  se = list(se1_iv_clust, se2_iv_clust, se3_iv_clust),
  t.auto = TRUE,
  no.space=TRUE,
  covariate.labels=c("Civilians killed by Govt, \\% of civilian deaths", "CPN(M) vote share",  
    "Madhesi Parties Vote Share", "Total Fatalities", 
    "Change in Economic Activity 1996-2006",
    "Economic activity 2006", "Average Slope of Area", 
    "Average Elevation of Area",
    "Constant"),
  omit.stat=c("ser", "f"),
  title="2SLS: Per capita aid by municipality (log USD)",
  label="tab:main_results_iv",
  table.placement = "!htb",
  float = FALSE,
#  notes="Robust standard errors. *p$<$0.1; **p$<$0.05; ***p$<$0.001.",
#  notes.append=FALSE,
#  notes.align = "r"
  add.lines = list(c(paste(rownames(summ_tot$diagnostics)[1], ", F-statistic", sep=""), 
                             signif(summ_tot$diagnostics[1, "statistic"], 2), 
                             signif(summ_soc$diagnostics[1, "statistic"], 2),
                             signif(summ_infr$diagnostics[1, "statistic"], 2)), 
                           c(paste(rownames(summ_tot$diagnostics)[2], ", p-value", sep=""), 
                             signif(summ_tot$diagnostics[2, "p-value"], 2), 
                             signif(summ_soc$diagnostics[2, "p-value"], 2),
                             signif(summ_infr$diagnostics[2, "p-value"],2)))
)

star_final_iv_clust <- star_insert_row(star_iv_clust, "\\textbf{\\textit{Political Support}}\\\\", insert.after = 11)
star_final_iv_clust <- star_insert_row(star_final_iv_clust, "\\textbf{\\textit{Controls: Reconstruction Need}}\\\\", insert.after = 18)
star_final_iv_clust <- star_insert_row(star_final_iv_clust, "\\textbf{\\textit{Controls: Economic Need}}\\\\", insert.after = 23)


cat(star_final_iv_clust, file = "main_aid_iv_clust.tex", sep="\n")


#########
#Appendix 

#########
#Table A1 -- Bypassing

#prep data
allaid_nep <- crs_raw %>% 
  filter(RecipientName == "Nepal", 
    Year >= 2006 & Year <= 2013) %>%
  group_by(Year) %>%
  summarize(aid_commited = sum(USD_Commitment_Defl, na.rm=TRUE),
    aid_disbursed = sum(USD_Disbursement_Defl, na.rm=TRUE))

govchan_nep <- crs_raw %>%
  filter(RecipientName == "Nepal", 
    Year >= 2006 & Year <= 2013) %>%
  # filter(ParentChannelCode %in% c(10000, 11000, 12000, 13000)) %>%
  filter(ParentChannelCode %in% c(10000, 12000, 13000)) %>%
  group_by(Year) %>%
  summarize(aid_commited_gov = sum(USD_Commitment_Defl, na.rm=TRUE),
    aid_disbursed_gov = sum(USD_Disbursement_Defl, na.rm=TRUE))

crs_nep <- allaid_nep %>%
  full_join(govchan_nep) %>%
  mutate(byp_com = 1-(aid_commited_gov/aid_commited),
    byp_disp = 1-(aid_disbursed_gov/aid_disbursed))

#create table
bypas_mat <- matrix(c(crs_nep$Year, crs_nep$byp_com),
  8, 2)
	#add row and column names
colnames(bypas_mat) <- c("Year", "Bypass, %")

tab_bypas <- stargazer(bypas_mat,
  style = "ajps",
  table.placement = "!htb",
  float = FALSE)

cat(tab_bypas, file = "table_bypas.tex", sep="\n")


#########
#Tables A2 & A3-- Variance inflation factors,
#Models with UML included

#Estimate all models 
##### All aid
F_TotAid_ln_elect_UML_iv <- formula(aid_all_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +

  Vote.share.UML +

  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean |
  CivGovShare | N_battalions_50)


M_TotAid_ln_elect_UML_iv <- ivreg(F_TotAid_ln_elect_UML_iv, data=df_nepal) 

## Social aid
F_SocAid_ln_elect_UML_iv <- formula(aid_socl_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +

  Vote.share.UML +

  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean |
  CivGovShare | N_battalions_50)


M_SocAid_ln_elect_UML_iv <- ivreg(F_SocAid_ln_elect_UML_iv, data=df_nepal) 

### Infrastructure aid
F_InfrAid_ln_elect_UML_iv <- formula(aid_infr_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +

  Vote.share.UML +

  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean |
  CivGovShare | N_battalions_50)

M_InfrAid_ln_elect_UML_iv <- ivreg(F_InfrAid_ln_elect_UML_iv, data=df_nepal) 

se1_UML_iv <- sqrt(diag(vcovCL(M_TotAid_ln_elect_UML_iv, cluster= ~NAME_2, type = "HC1")))
se2_UML_iv <- sqrt(diag(vcovCL(M_SocAid_ln_elect_UML_iv, cluster= ~NAME_2, type = "HC1")))
se3_UML_iv <- sqrt(diag(vcovCL(M_InfrAid_ln_elect_UML_iv, cluster= ~NAME_2, type = "HC1")))

#Table A3: Variance Inflation Facto; Note: make sure to also run main models first!
vif_mat <- matrix(c(vif(M_TotAid_ln_elect_iv)[2:3], NA,
  vif(M_TotAid_ln_elect_UML_iv)[2:4]), 3, 2)

#Now add row and column names
rownames(vif_mat) <- c("Maoist Vote share", "Madhesi Parties Vote Share",
  "UML Vote Share")
colnames(vif_mat) <- c("Restricted", "UML included")

tab_vif <- stargazer(vif_mat,
  title = "Variance inflation factors: IV model, restricted and with UML vote share included.",
  label = "tab:vif",
  style = "ajps",
  table.placement = "!htb",
  float = FALSE)

cat(tab_vif, file = "table_vif.tex", sep="\n")


#########
#Table A4

#manualy remove observations with missing values, using spreg results
wfull <- nb2mat(nb, style="B", zero.policy=TRUE)
wmod1 <- wfull[-M_TotAid_SAC$na.action,-M_TotAid_SAC$na.action]
wmod2 <- wfull[-M_SocAid_SAC$na.action, -M_SocAid_SAC$na.action]
wmod3 <- wfull[-M_InfrAid_SAC$na.action, -M_InfrAid_SAC$na.action]
lw1 <- mat2listw(wmod1, style="W")
lw2 <- mat2listw(wmod2, style="W")
lw3 <- mat2listw(wmod3, style="W")


table_main_UML_iv <- list(M_TotAid_ln_elect_UML_iv, M_SocAid_ln_elect_UML_iv, 
  M_InfrAid_ln_elect_UML_iv)


  #collect diagnostic statistics
summ_tot <- summary(M_TotAid_ln_elect_UML_iv, vcov. = function(x) 
  vcovCL(x, cluster =df_nepal[-M_TotAid_ln_elect_UML_iv$na.action,]$NAME_2, type="HC1"), diagnostics=T)
summ_soc <- summary(M_SocAid_ln_elect_UML_iv, vcov. = function(x) 
  vcovCL(x, cluster = df_nepal[-M_SocAid_ln_elect_UML_iv$na.action,]$NAME_2, type="HC1"), diagnostics=T)
summ_infr <- summary(M_InfrAid_ln_elect_UML_iv, vcov. = function(x) 
  vcovCL(x, cluster = df_nepal[-M_InfrAid_ln_elect_UML_iv$na.action,]$NAME_2, type="HC1"), diagnostics=T)

star_main_UML_iv <- stargazer(table_main_UML_iv,
  # order = c(1, 7, 8, 2, ),
  dep.var.labels=c("Total Aid", "Social Aid", "Infrastructure\\\\&&&Aid"),
  se = list(se1_UML_iv, se2_UML_iv, se3_UML_iv),
  t.auto = TRUE,
  no.space=TRUE,
  covariate.labels=c("Civilians killed by Govt, \\% of civilian deaths", "CPN(M) vote share",  
    "Madhesi Parties Vote Share", 
    "UML Vote Share", 
  "Total Fatalities", 
    "Change in Economic Activity 1996-2006",
    "Economic activity 2006", "Average Slope of Area", 
    "Average Elevation of Area",
    "Constant"),
  omit.stat=c("ser", "f"),
  title="2SLS: Per capita aid by municipality (log USD)",
  label="tab:main_results_UML_iv",
  table.placement = "!htb",
  float = FALSE,
#  notes="Robust standard errors. *p$<$0.1; **p$<$0.05; ***p$<$0.001.",
#  notes.append=FALSE,
#  notes.align = "r"
  add.lines = list(c(paste(rownames(summ_tot$diagnostics)[1], ", F-statistic", sep=""), 
                             signif(summ_tot$diagnostics[1, "statistic"], 2), 
                             signif(summ_soc$diagnostics[1, "statistic"], 2),
                             signif(summ_infr$diagnostics[1, "statistic"], 2)), 
                           c(paste(rownames(summ_tot$diagnostics)[2], ", p-value", sep=""), 
                             signif(summ_tot$diagnostics[2, "p-value"], 2), 
                             signif(summ_soc$diagnostics[2, "p-value"], 2),
                             signif(summ_infr$diagnostics[2, "p-value"],2)))
)

star_main_final_UML_iv <- star_insert_row(star_main_UML_iv, "\\textbf{\\textit{Political Support}}\\\\", insert.after = 11)
star_main_final_UML_iv <- star_insert_row(star_main_final_UML_iv, "\\textbf{\\textit{Controls: Reconstruction Need}}\\\\", insert.after = 20)
star_main_final_UML_iv <- star_insert_row(star_main_final_UML_iv, "\\textbf{\\textit{Controls: Economic Need}}\\\\", insert.after = 25)

cat(star_main_final_UML_iv, file = "main_aid_iv_UML.tex", sep="\n")


#########
#Table A5

#declare panel data 
pdata <- pdata.frame(df_nepal, index = "NAME_2")

#estimate models

#total aid, IV
M_TotAid_IV_plm <- plm(aid_all_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean +
  CivGovShare | N_battalions_50 + 
   Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean ,
  data = pdata,
  model = "random",
  random.method = "walhus")

#social aid, IV
M_SocAid_IV_plm <- plm(aid_socl_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean +
  CivGovShare | N_battalions_50 + 
   Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean ,
  data = pdata,
  model = "random",
  random.method = "walhus")

#infrastructure aid, IV
M_InfrAid_IV_plm <- plm(aid_infr_tot_pc_ln ~
  #Political variables
  Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean +
  CivGovShare | N_battalions_50 + 
   Vote.share.Maoist + I(Vote.share.Terai_madesh+Vote.share.Madeshi_jana) +
  #Reconstruction needs
  totkill + I(v4composites_calibrated_201709.2006.mean - v4composites_calibrated_201709.1996.mean) +
  #Economic need & fighting difficulty
  v4composites_calibrated_201709.2006.mean + srtm_slope_500m.none.mean + srtm_elevation_500m.none.mean ,
  data = pdata,
  model = "random",
  random.method = "walhus")

se1_re_IV_clust <- sqrt(diag(plm::vcovHC(M_TotAid_IV_plm, cluster= "group", type = "HC1", method = "white2")))
se2_re_IV_clust <- sqrt(diag(plm::vcovHC(M_SocAid_IV_plm, cluster= "group", type = "HC1", method = "white2")))
se3_re_IV_clust <- sqrt(diag(plm::vcovHC(M_InfrAid_IV_plm, cluster= "group", type = "HC1", method = "white2")))


#Create table
table_re_iv <- list(M_TotAid_IV_plm, M_SocAid_IV_plm, 
  M_InfrAid_IV_plm)

star_re_iv <- stargazer(table_re_iv,
  dep.var.labels=c("Total Aid", "Social Aid", "Infrastructure\\\\&&&Aid"),
  se = list(se1_re_IV_clust, se2_re_IV_clust, se3_re_IV_clust),
  t.auto = TRUE,
  no.space=TRUE,
  order =c(8, 1,2,3,4,5,6,7,9),
  covariate.labels=c(    "Civilians killed by Govt, \\% of civilian deaths",
    "CPN(M) vote share",  
    "Madhesi Parties Vote Share", "Total Fatalities", 
    "Change in Economic Activity 1996-2006",
    "Economic activity 2006", "Average Slope of Area", 
    "Average Elevation of Area",
    "Constant"),
  omit.stat=c("ser", "f"),
  title="OLS with district-random effects: Per capita aid by municipality (log USD)",
  label="tab:results_re",
  table.placement = "!htb",
  float = FALSE
#  notes="Robust standard errors. *p$<$0.1; **p$<$0.05; ***p$<$0.001.",
#  notes.append=FALSE,
#  notes.align = "r"
)

star_re_iv_final <- star_insert_row(star_re_iv, 
  "\\textbf{\\textit{Political Support}}\\\\", insert.after = 11)
star_re_iv_final <- star_insert_row(star_re_iv_final, 
  "\\textbf{\\textit{Controls: Reconstruction Need}}\\\\", insert.after = 18)
star_re_iv_final <- star_insert_row(star_re_iv_final, 
  "\\textbf{\\textit{Controls: Economic Need}}\\\\", insert.after = 23)

cat(star_re_iv_final, file = "main_re_iv.tex", sep="\n")


